In-class Exercise 5

A short description of the post.

Lye Jia Wei https://lye-jia-wei.github.io/
09-13-2021

Important to set fig.retina to be 3 in order for subsequent visualization to be sharp.

Installing and Loading the R package

packages <- c('maptools','sf','raster','spatstat','tmap','tidyverse', 'plotly','ggthemes')
for (p in packages){
  if(!require(p, character.only = T)){
    install.packages(p)
  }
}
library(p, character.only=T)

Import the Geospatial Data

The code chunk import shapefile using st_read() of sf package. The output object is in tibble sf object class.

Pgl_sf and mpsz_sf geospatial data sets into RStudio.

Note: - Do not paste directory in dsn

Pgl_sf <- st_read(dsn="data/shapefile", layer="Punggol_St")
Reading layer `Punggol_St' from data source 
  `C:\Users\User\Desktop\IS415\lye-jia-wei\IS415_blog\_posts\2021-09-13-in-class-exercise-5\data\shapefile' 
  using driver `ESRI Shapefile'
Simple feature collection with 2642 features and 2 fields
Geometry type: LINESTRING
Dimension:     XY
Bounding box:  xmin: 34038.56 ymin: 40941.11 xmax: 38882.85 ymax: 44801.27
Projected CRS: SVY21 / Singapore TM
mpsz_sf <- st_read(dsn="data/shapefile",layer="MP14_SUBZONE_WEB_PL")
Reading layer `MP14_SUBZONE_WEB_PL' from data source 
  `C:\Users\User\Desktop\IS415\lye-jia-wei\IS415_blog\_posts\2021-09-13-in-class-exercise-5\data\shapefile' 
  using driver `ESRI Shapefile'
Simple feature collection with 323 features and 15 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 2667.538 ymin: 15748.72 xmax: 56396.44 ymax: 50256.33
Projected CRS: SVY21

Projection is SVY21

Import aspatial data from rds folder

read_rds() of readr package is used instead of readRDS() of base R is used. This is because output of read_rds() is in tibble object.

Note:

childcare <- read_rds("data/rds/childcare.rds")
CHAS <- read_rds("data/rds/CHAS.rds")

Converting the aspatial data frame into sf objects

CHAS_sf <- st_as_sf(CHAS, coords = c("X_COORDINATE","Y_COORDINATE"),crs=3414)


childcare_sf <- st_as_sf(childcare, coords= c("Lng","Lat"),crs=4326) %>% st_transform(crs=3414)

Note:

Plotting for Reviewing

Alpha is to set transparency

tmap_mode("view")
tm_shape(childcare_sf)+tm_dots(alpha= 0.4, col="blue", size=0.05)+ tm_shape(CHAS_sf)+tm_dots(alpha= 0.4, col="red", size=0.05)

Note:

Geospatial Data Wrangling

Converting from sf to Spatial* data frame

as_Spatial() of sf package.

childcare <- as_Spatial(childcare_sf)
CHAS <- as_Spatial(CHAS_sf)
mpsz <- as_Spatial(mpsz_sf)

Converting Spatial* data frame into Spatial* Objects

as.SpatialPoint() of as.SpatialPolygon() of maptools package

childcare_sp <- as(childcare, "SpatialPoints")
CHAS_sp <- as(CHAS, "SpatialPoints")
MPSZ_sp <- as(mpsz, "SpatialPolygons")

Converting from Spatial* objects into ppp objects

This is to drop the geographical projection information using as.ppp() of maptools package.

childcare_ppp <- as(childcare_sp, "ppp")
CHAS_ppp <- as(CHAS_sp, "ppp")

Removing duplicate points using jitter

childcare_ppp_jit <- rjitter(childcare_ppp, retry = TRUE, nsim=1, drop = TRUE)

any(duplicated(childcare_ppp_jit))  
[1] FALSE

Output return a ‘False’ which suggest that there are no more duplicate.

Note:

CHAS_ppp_jit <- rjitter(CHAS_ppp, retry = TRUE, nsim=1, drop = TRUE)

any(duplicated(CHAS_ppp_jit))  
[1] FALSE

Output return a ‘False’ which suggest that there are no more duplicate.

Note:

Extracting Punggol Planning Area

Extracting planning sub zone PUNGGOL from data and then planning area.

pg <- mpsz[mpsz@data$PLN_AREA_N=="PUNGGOL",]

Note:

Converting Spatial Polygon Data frame into Spatial Polygon object

pg_sp <- as(pg, "SpatialPolygons")

Converting Spatial Polygons into owin object

Convert into owin object which is a list of length 5.

pg_owin <- as(pg_sp, "owin")

Extracting spatial points window owin

childcare_pg <- childcare_ppp_jit[pg_owin]
CHAS_pg <- CHAS_ppp_jit[pg_owin]
plot(childcare_pg)

L-function

L_childcare <- envelope(childcare_pg,
                        Lest,
                        nsim = 99,
                        rank = 1,
                        gloval = TRUE)
Generating 99 simulations of CSR  ...
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35,
36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70,
71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98,  99.

Done.
title <- "Pairwise Distance: L function"

Lcsr_df <- as.data.frame(L_childcare)

colour=c("#0D657D","#ee770d","#D3D3D3")
csr_plot <- ggplot(Lcsr_df, aes(r, obs-r))+
  # plot observed value
  geom_line(colour=c("#4d4d4d"))+
  geom_line(aes(r,theo-r), colour="red", linetype = "dashed")+
  # plot simulation envelopes
  geom_ribbon(aes(ymin=lo-r,ymax=hi-r),alpha=0.1, colour=c("#91bfdb")) +
  xlab("Distance r (m)") +
  ylab("L(r)-r") +
  geom_rug(data=Lcsr_df[Lcsr_df$obs > Lcsr_df$hi,], sides="b", colour=colour[1])  +
  geom_rug(data=Lcsr_df[Lcsr_df$obs < Lcsr_df$lo,], sides="b", colour=colour[2]) +
  geom_rug(data=Lcsr_df[Lcsr_df$obs >= Lcsr_df$lo & Lcsr_df$obs <= Lcsr_df$hi,], sides="b", color=colour[3]) +
  theme_tufte()+
  ggtitle(title)

text1<-"Significant clustering"
text2<-"Significant segregation"
text3<-"Not significant clustering/segregation"

# the below conditional statement is required to ensure that the labels (text1/2/3) are assigned to the correct traces
if (nrow(Lcsr_df[Lcsr_df$obs > Lcsr_df$hi,])==0){ 
  if (nrow(Lcsr_df[Lcsr_df$obs < Lcsr_df$lo,])==0){ 
    ggplotly(csr_plot, dynamicTicks=T) %>%
      style(text = text3, traces = 4) %>%
      rangeslider() 
  }else if (nrow(Lcsr_df[Lcsr_df$obs >= Lcsr_df$lo & Lcsr_df$obs <= Lcsr_df$hi,])==0){ 
    ggplotly(csr_plot, dynamicTicks=T) %>%
      style(text = text2, traces = 4) %>%
      rangeslider() 
  }else {
    ggplotly(csr_plot, dynamicTicks=T) %>%
      style(text = text2, traces = 4) %>%
      style(text = text3, traces = 5) %>%
      rangeslider() 
  }
} else if (nrow(Lcsr_df[Lcsr_df$obs < Lcsr_df$lo,])==0){
  if (nrow(Lcsr_df[Lcsr_df$obs >= Lcsr_df$lo & Lcsr_df$obs <= Lcsr_df$hi,])==0){
    ggplotly(csr_plot, dynamicTicks=T) %>%
      style(text = text1, traces = 4) %>%
      rangeslider() 
  } else{
    ggplotly(csr_plot, dynamicTicks=T) %>%
      style(text = text1, traces = 4) %>%
      style(text = text3, traces = 5) %>%
      rangeslider()
  }
} else{
  ggplotly(csr_plot, dynamicTicks=T) %>%
    style(text = text1, traces = 4) %>%
    style(text = text2, traces = 5) %>%
    style(text = text3, traces = 6) %>%
    rangeslider()
  }